Load network and scent data

From Tables S1 and S2 in: Kantsa, A., R.A. Raguso, A.G. Dyer, J.M. Olesen, T. Tscheulin, and T. Petanidou. 2018. Disentangling the role of floral sensory stimuli in pollination networks. Nature Communications 9: 1041. doi:[10.1038/s41467-018-03448-w](https://doi.org/10.1038/s41467-018-03448-w)

#load network
net <- read.delim("Kantsa_etal_Supp_Data1.csv", skip = 2)
rownames(net) <- net[,1]
net[,1] <- NULL
net <- as.matrix(t(net))
snet <- sortweb(net) #sort by row & column totals

#load floral scents
scent.file <- "Kantsa_2018_Supp_Data2.xlsx"
plants <- getSheetNames(scent.file)
scentlist <- lapply(2:length(plants), read.xlsx, xlsxFile=scent.file)
plants <- plants[-1]
scent <- rbindlist(scentlist, fill=T, idcol="plant")[,-4]
colnames(scent)[2:3] <- c("cmpd","mean")
scent$plant <- factor(plants[scent$plant])
scent$cmpd <- factor(tolower(scent$cmpd))
cmpds <- levels(scent$cmpd)

scent.net <- dcast(scent, plant ~ cmpd, value.var="mean", fun.aggregate=mean, fill=0) #long to wide
scent.net <- as.data.frame(scent.net)
rownames(scent.net) <- scent.net[,1]
scent.net[,c(1,380)] <- NULL #weird NA column has one entry
scent.net <- as.matrix(scent.net)

Plant phylogeny

#plants.class <- classification(plants, db="gbif")
#save(plants.class, file="plantsclass.Rdata")
load("plantsclass.Rdata")
plants.tree <- class2tree(plants.class)
plot(plants.tree)

Plant-pollinator network

Heatmap

heatmaply(net, scale="column")

Matrix plot

#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

Web plot

#interaction web plot
plotweb(net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])

Network topography

networklevel(net)
                         connectance                     web asymmetry 
                          0.06312657                        0.63106796 
                   links per species            number of compartments 
                          1.95631068                        2.00000000 
               compartment diversity               cluster coefficient 
                          1.07901412                        0.02631579 
                          nestedness               weighted nestedness 
                          4.48059576                        0.58672881 
                       weighted NODF    interaction strength asymmetry 
                         10.80348739                        0.20372363 
            specialisation asymmetry                   linkage density 
                         -0.04424002                        3.57255189 
                weighted connectance                      Fisher alpha 
                          0.01734248                       89.61252715 
                   Shannon diversity              interaction evenness 
                          3.27629195                        0.37393976 
        Alatalo interaction evenness                                H2 
                          0.25285628                        0.57919894 
                number.of.species.HL              number.of.species.LL 
                        168.00000000                       38.00000000 
   mean.number.of.shared.partners.HL mean.number.of.shared.partners.LL 
                          0.34944397                        1.18776671 
              cluster.coefficient.HL            cluster.coefficient.LL 
                          0.34062967                        0.21733898 
     weighted.cluster.coefficient.HL   weighted.cluster.coefficient.LL 
                          0.72421977                        0.31277076 
                    niche.overlap.HL                  niche.overlap.LL 
                          0.12366249                        0.14889538 
                     togetherness.HL                   togetherness.LL 
                          0.03019814                        0.02302573 
                          C.score.HL                        C.score.LL 
                          0.75923950                        0.74100046 
                          V.ratio.HL                        V.ratio.LL 
                          3.23242952                       17.02718828 
                      discrepancy.HL                    discrepancy.LL 
                        248.00000000                      250.00000000 
                 extinction.slope.HL               extinction.slope.LL 
                          5.03260774                        1.60178636 
                       robustness.HL                     robustness.LL 
                          0.81123903                        0.61539429 
       functional.complementarity.HL     functional.complementarity.LL 
                       8044.31456516                     7651.23911634 
                partner.diversity.HL              partner.diversity.LL 
                          0.94103312                        1.28483306 
                       generality.HL                  vulnerability.LL 
                          2.67434342                        4.47076037

Scent network

Heatmap

heatmaply(scent.net, scale="column")

Web plot

#interaction web plot
plotweb(scent.net, text.rot=90, method="cca", col.high=pal[2], col.low=pal[1], col.interaction=pal[5])

NMDS

scent.mds <- metaMDS(sqrt(scent.net), autotransform = F)
   Run 0 stress 0.2177614 
   Run 1 stress 0.2155518 
   ... New best solution
   ... Procrustes: rmse 0.08183959  max resid 0.3861441 
   Run 2 stress 0.2277 
   Run 3 stress 0.2198936 
   Run 4 stress 0.2217743 
   Run 5 stress 0.2139972 
   ... New best solution
   ... Procrustes: rmse 0.1078209  max resid 0.3550252 
   Run 6 stress 0.2201085 
   Run 7 stress 0.2306285 
   Run 8 stress 0.2158631 
   Run 9 stress 0.2417021 
   Run 10 stress 0.2264037 
   Run 11 stress 0.2330208 
   Run 12 stress 0.2175837 
   Run 13 stress 0.2506198 
   Run 14 stress 0.2138077 
   ... New best solution
   ... Procrustes: rmse 0.07751348  max resid 0.2765495 
   Run 15 stress 0.2130163 
   ... New best solution
   ... Procrustes: rmse 0.08566332  max resid 0.3435179 
   Run 16 stress 0.2426195 
   Run 17 stress 0.2136397 
   Run 18 stress 0.2156501 
   Run 19 stress 0.2146277 
   Run 20 stress 0.2169694 
   *** No convergence -- monoMDS stopping criteria:
       20: stress ratio > sratmax
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col=pal[3])
text(scent.op, what="sites", cex=0.8)

Phylochemospace

scent.mds.taxa <- scent.mds$points
rownames(scent.mds.taxa) <- plants.tree$phylo$tip.label

plants.tree$phylo$edge.length <- plants.tree$phylo$edge.length+0.00001

phylomorphospace(plants.tree$phylo, scent.mds.taxa, control=list(col.node=0), label="none", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T)

Phylo Bar chart

par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(scent.net), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE) 

par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1)
par(mar=c(0,0,0,0))
barplot(t(decostand(scent.net, method="tot")), col=sample(rainbow(ncol(scent.net))), names.arg=rep("",nrow(scent.net)), horiz=TRUE)